rk = [];
for alpha=0.01:0.01:0.99
    N=100;
    h=0.01;
    z=sqrt(3)/pi*(log(alpha)-log(1-alpha));
    y(1)=h;
    for i=1:N
        K1=y(i)*(1+z)+1;
        K2=y(i)*(1+z)+1+h/2*(y(i)*(1+z)+1)*(1+z);
        K3=y(i)*(1+z)+1+h/2*(y(i)*(1+z)+1)*(1+z)+h/4*(y(i)*(1+z)+1)*(1+z)*(1+z);
        K4=y(i)*(1+z)+1+h*(y(i)*(1+z)+1)*(1+z)+h^2/2*(y(i)*(1+z)+1)*(1+z)*(1+z)+h^2/4*(y(i)*(1+z)+1)*(1+z)*(1+z)*(1+z);
        y(i+1)=y(i)+h/6*(K1+2*K2+2*K3+K4);
    end
    rk = [rk;y(i)];
end
plot(rk,0.01:0.01:0.99,'g-.')

